//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : timederivatives_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef TIMEDERIVATIVES_KERNELS_H
#define TIMEDERIVATIVES_KERNELS_H

#ifdef USEDOUBLE
#define ONE 1.0
#define HALF 0.5
#else
#define ONE 1.0f
#define HALF 0.5f
#endif

namespace wavescomputer
{

template <class T>
__global__ void kernel_timederivatives(T *dt_zeta,
                                       T *dt_phi,
                                       const T *U1,
                                       const T *U2,
                                       const T *h0,
                                       const T *D0,
                                       const T *zeta,
                                       const T *phi,
                                       const T *psi,
                                       const T *SRCz,
                                       const T *pressurePulse,
                                       const T pslow,
                                       const int Nx1,
                                       const int Nx2,
                                       const T dx1,    // mesh spacing x
                                       const T dx2,    // mesh spacing y
                                       const int ld,
                                       const T g)
{
    int u = BORDER_WIDTH + bx * Bx + tx;
    int v = BORDER_WIDTH + by * By + ty;
    int tid = v * ld + u;  

    const T dx1invtwice = (HALF / dx1);
    const T dx2invtwice = (HALF / dx2);
    const T dx1invsquared = (ONE / (dx1 * dx1));
    const T dx2invsquared = (ONE / (dx2 * dx2));

    const T U1C = U1[tid];
    const T U2C = U2[tid];
    const T SzzE = - dx1invtwice * (HALF * (U1C + U1[tid + ld]));
    const T SzzW =   dx1invtwice * (HALF * (U1C + U1[tid - ld]));
    const T SzzN = - dx2invtwice * (HALF * (U2C + U2[tid + 1]));
    const T SzzS =   dx2invtwice * (HALF * (U2C + U2[tid - 1]));
    const T SzzC = SzzE + SzzW + SzzN + SzzS;

    const T h0C = h0[tid];
    const T h0Ebar = HALF * (h0C + h0[tid + ld]);
    const T h0Wbar = HALF * (h0C + h0[tid - ld]);
    const T h0Nbar = HALF * (h0C + h0[tid + 1]);
    const T h0Sbar = HALF * (h0C + h0[tid - 1]);

    const T SzfE = dx1invsquared * h0Ebar;
    const T SzfW = dx1invsquared * h0Wbar;
    const T SzfN = dx2invsquared * h0Nbar;
    const T SzfS = dx2invsquared * h0Sbar;
    const T SzfC = SzfE + SzfW + SzfN + SzfS;

    const T D0C = D0[tid];
    const T SzpE = dx1invsquared * h0Ebar * (HALF * (D0C + D0[tid + ld]));
    const T SzpW = dx1invsquared * h0Wbar * (HALF * (D0C + D0[tid - ld]));
    const T SzpN = dx2invsquared * h0Nbar * (HALF * (D0C + D0[tid + 1]));
    const T SzpS = dx2invsquared * h0Sbar * (HALF * (D0C + D0[tid - 1]));
    const T SzpC = SzpE + SzpW + SzpN + SzpS;

    const T phiE = phi[tid + ld];
    const T phiW = phi[tid - ld];
    const T phiN = phi[tid + 1 ];
    const T phiS = phi[tid - 1 ];
    const T phiC = phi[tid];
    const T zetaC = zeta[tid];

    dt_zeta[tid] =   SzzE * zeta[tid + ld]
                   + SzzW * zeta[tid - ld]
                   + SzzN * zeta[tid + 1]
                   + SzzS * zeta[tid - 1]
                   - SzzC * zetaC
                   - SzfE * phiE
                   - SzfW * phiW
                   - SzfN * phiN
                   - SzfS * phiS
                   + SzfC * phiC
                   + SzpE * psi[tid + ld]
                   + SzpW * psi[tid - ld]
                   + SzpN * psi[tid + 1]
                   + SzpS * psi[tid - 1]
                   - SzpC * psi[tid]
                   - pslow * SRCz[tid];

     dt_phi[tid] = - g * zetaC
                   + SzzE * phiE
                   + SzzW * phiW
                   + SzzN * phiN
                   + SzzS * phiS
                   - SzzC * phiC
                   - pressurePulse[tid];
}

} // namespace wavescomputer

#endif // TIMEDERIVATIVES_KERNELS_H

